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Abstract 

The phase behavior of colloidal particles embedded in a binary fluid is in- 
fluenced by wetting layers surrounding each particle. The free energy of the 
fluid film depends on its morphology, i.e., on size, shape and connectivity. Un- 
der rather natural assumptions a general expression for the Hamiltonian can 
be given extending the model of hard spheres to partially penetrable shapes 
including energy contributions related to the volume, surface area, mean cur- 
vature, and Euler characteristic of the wetting layer. The complex spatial 
structure leads to multi-particle interactions of the colloidal particles. 

The dependence of the morphology of the wetting layer on temperature and 
density can be studied using Monte-Carlo simulations and perturbation theory. 
A fluid-fluid phase separation induced by the wetting layer is observed which is 
suppressed when only two-particle interactions are taken into account instead 
of the inherent many-particle interaction of the wetting layer. 

Key words: colloids, hard-disks, integral geometry, Monte-Carlo simulation, 
penetrable spheres, phase-transition, topology, wetting. 



1 Introduction 



Colloidal solutions such as paints and soots are extremely common and exhibit many 
industrial and technological applications. Due to their mesoscopic size they are also 
of great fundamental importance and have provided basic parts of our understanding 
of the interaction of particles, for instance. In most cases the electrostatic repulsive 
interaction is significantly screened so that the attractive dispersion interaction dom- 
inate which finally cause an irreversible aggregation or coagulation. The colloidal 
particles stick, i.e., are in point contact at the global minimum of the interaction 
potential. But under certain circumstance reversible colloidal aggregation or floc- 
culation can be observed. Systematic experiments [|l], 0, |3| on colloidal particles 
embedded in a near-critical solvent mixture of 2,6-lutidine plus water, for instance, 
have revealed that flocculation can be viewed as thermally induced phase separation. 

Theoretical work of colloidal partitioning or flocculation in a two-phase solvent fo- 
cused on general thermodynamic arguments known from wetting or capillary conden- 
sation [|, ^| and more recently on a Ginzburg-Landau approach including a detailed 
treatment of the first-order solvent phase transition which is coupled to the colloidal 
particles || 0. These works are based on microscopic pair interactions between the 
colloidal particles and the molecules in the binary fluid. Thus, the colloids do not 
interact directly, but couple to the binary fluid degrees of freedom by preferring one 
of the two components. 

The phase behavior of colloidal particles embedded in a binary fluid is certainly 
influenced by the existence of wetting layers, i.e. by thermodynamically meta-stable 
fluid phases stabilized at the boundary of the colloidal particles. Wetting phenomena 
appear in multicomponent systems when the components exhibit different interac- 
tions with the colloidal particle. In general, when two thermodynamic phases are 
close to a first order phase transition the wetting layer may become large and compa- 
rable to the diameter of the colloidal particle. The interactions are then determined 



by the free energy of the fluid film between the hard colloidal particles which cause 
clustering and eventually a phase separation. 

Here, we focus on an effective theory where the binary fluid enters only in one 
parameter, namely the thickness L(T, c) of the wetting layer (enriched layer of one 
component) which depends on temperature T and concentration c of one component 
of the binary fluid. In other words, the colloidal system and the binary fluid are cou- 
pled solely by the parameter L(T, c) which can be determined for a single colloidal 
particle in a binary fluid without taking into account the interaction with other parti- 
cles. The advantage of this approach is the computational simplicity of the influence 
of the binary fluid on the interaction of the colloidal particles. This simplicity allows 
a more sophisticated treatment of the colloidal interaction itself. Taking into account 
not only pair-potentials one can study, for instance, the effect of multiple-particle 
interactions due to an overlap of the wetting layer of many particles. This cannot be 
neglected when L is large in the complete wetting regime, i.e., close to the critical 
point of the binary fluid. 

Since the thermodynamic properties of the fluid film depend on its morphology, 
i.e., on volume and surface statistical theory should include geometrical de- 

scriptors to characterize the size, shape and connectivity of the wetting layer. In 
Section || we proposed a model for the study of such colloidal suspensions. The col- 
loids are resembled by spherical particles (disks in two dimensions) with a hard-core 
diameter D and a soft, penetrable shell of thickness L (see Figure [I]). Using the radius 
R = D/2 + L of the partially penetrable disks one can define the ratio 5 = D/(2R) 
with < 5 < 1, where 5 = 1 equals a pure hard-disk system, whereas 5 = denotes 
an ensemble of fully penetrable disks. The Hamiltonian includes not only the two- 
particle hard-sphere potential but also energy contributions related to the volume, 
surface area, mean curvature, and Euler characteristic of the layers around the hard 
particles. 

In Section |] a perturbation theory and in Section f| a Metropolis Monte-Carlo al- 



gorithm of partially penetrable disks in two dimensions is presented. The morphology 
of the wetting layer as function of temperature and density of the colloidal particles 
as well as phase diagrams are discussed. 

2 Morphological thermodynamics of colloidal con- 
figurations 

From a geometric point of view one should tie the interaction of the colloidal particles 
to the morphology of the wetting layer, for instance, to its volume and surface area. 
Each configuration 

N 

^ = U^) • (i) 

1=1 

is assumed to be the union of mutually penetrable ci-dimensional spheres B^(x) of 
radius R and hard-core diameter D centered at x G O C embedded in binary 
fluid host component. For convenience, we assume a box Q of edge length 1 and 
periodic boundary conditions so that all lengths are measured in units of the box 
size. Typical configurations are shown in Figure [l] for 5 = D/(2R) = 0.4. In the 
following we consider two dimensional systems with penetrable disks. Depending 
on the density p of the particles and the size L of the soft shell, i.e. the ratio 
5 = 1/(1 + 2L/D) of the radii, the wetting layer (i.e., the white area in Figure 
exhibit quit different topological and geometric properties. For instance, the white 
disk shells are disconnected (isolated) for S — > 1 due to the hard core interaction, 
whereas for small S ~ the grains can overlap so that at higher densities p connected 
structures occur. The morphology of the emerging pattern may be characterized by 
the covered volume or the area in two dimensions, the surface area or boundary length, 
respectively, and the Euler-characteristic, i.e., the connectivity of the penetrating 
grains. Area, boundary length and Euler characteristic have in common that they 
are additive measures. Additivity means that the measure of the union A U B of 



two domains (grains) A and B equals the sum of the measure of the single domains 
subtracted by the intersection, i.e., 

M(AUB) = M{A)+M(B)-M(AnB) . (2) 

A remarkable theorem in integral geometry || |9|] is the completeness of the so-called 
Minkowski functionals. The theorem asserts that any additive, motion invariant and 
conditional continuous functional A4(A) = ^2^, =0 (h/M v (A) on subsets A C M d , is a 
linear combination of the d + 1 Minkowski functionals M u , {i> — 0, . . . , d) with real 
coefficients c v independent of A. The two conditions of motion invariance and con- 
ditional continuity are necessary for the theorem, but they are not very restrictive in 
most physical situations. Intuitively, conditional continuity expresses the fact that 
an approximation of a convex domain K by convex polyhedra K n , for example, also 
yields an approximation of A4(K) by Ai(K n ). The property of motion-invariance 
means that the morphological measure M. of a domain A is independent of its lo- 
cation and orientation in space. Thus, every morphological measure M. which is 
additive (motion-invariant and continuous), i.e., which obeys relation (|2]) can be writ- 
ten in terms of Minkowski functionals M v , which are related to curvature integrals 
and do not only characterize the size but also shape (morphology) and connectivity 
(topology) of spatial patterns. In the three-dimensional Euclidean space the family 
of Minkowski functionals consists of the volume V = Mo, the surface area S = 8 Mi 
of the pattern, its integral mean curvature H = 2ir 2 M 2 , and integral of Gaussian cur- 
vature. In two dimension they are given by the area F = Mq, the length U = 2nMi 
and the curvature integral, i.e., the Euler characteristic x — ttM 2 of the boundary. 
In other words, the Minkowski functionals are the complete set of additive, morpho- 
logical measures. We assume, that the energy of a configuration is a morphological 
measure of the wetting layer, i.e., an additive functional of the fluid films surrounding 



each particle. Thus, the Boltzmann weights are specified by the Hamiltonian 

h(a n ) = v &^i) + \ M » (U^(^)) - NM »( B £v\ ( 3 ) 

i^j=l u=0 \ \i=l J J 

which is a linear combination of Minkowski functionals of the configurations and a 
pair potential 

{0 \xi — Xj\ > D 
( 4 ) 
oo \xi — Xj\ < D 

between two colloidal particles located at Xi and Xj. For convenience we assume a 
pure hard-core interaction. We emphasize that the Hamiltonian (Eq. @) constitutes 
the most general model for composite media assuming additivity of the free energy of 
the homogeneous wetting layer. The interactions between the colloidal particles are 
given by a bulk term (volume energy), a surface term (surface tension), and curvature 
terms (bending energies) of the wetting layer. 

We define the packing fraction r\ = jD 2 p = x5 2 and the normalized density 
x = puR 2 of the disks. The closest packing fraction is r/cp = or xcp = Vcp/8 2 , 
i.e., xcp(S = 0.4) w 5.67. Since we are not interested in the solid-liquid phase 
transition of hard disks but in a fluid-fluid transition induced by the wetting layer 
we focus on densities well below xcp- Depending on the ratio 5 it is possible to have 
multiple-overlapps of disks. In particular, two disks can overlapp iff 5 < 1, three 
iff 5 < V3/2 ~ 0.866, and four iff 5 < 1/V3 ~ 0.577. If 5 < 1/2 a proliferation of 
possible overlapps occur. In this limit a virial expansion of the Minkowski functionals 
in the density is not useful anymore and the interaction has manifestly many-body 
character. Therefore, we focus in this paper on 5 = 0.4 in comparision with S = 1 
and 5 = as limiting cases. 



3 Perturbation theory 

Because of the proliferation of multibody potentials, an exact evaluation of the par- 
tition function for the Hamiltonian (Q) appears to be unmanagable for d > 2. For 



the sake of comparison with Monte-Carlo simulations presented in the next section, 
we calculate here the thermodynamic properties and phase diagram of the model 
by simple first-order thermodynamic perturbation theory fl2 |. This approximation 



keeps the geometrical and topological aspects of the model intact. As a reference 
system, we use the hard-sphere fluid and solid. The free energy is estimated as 

(3F{T, p) = (3F HS (p) +(3<U> HS (p) (5) 

where F H s(p) is the free energy of the hard-sphere system at density p, and < U >hs 
(p) is the average value of the morphological part of the Hamiltonian (|3|), computed 
in the hard-sphere reference system. Thus, keeping only the first two terms in a high- 
temperature expansion of the free energy amounts to replacing the configurational 
integral in the partition function by exp {—(3 < U >hs} which yields a lower bound. 
The equation of state for hard disks is given by the Pade approximation from Hoover 



and Ree (1969) for the free energy [13| 



N B KH ' H 1 - 0.67bp + 0.096V 

where A is the mean thermal de Broglie wavelength and b = \D 2 . The packing 
fraction above the freezing transition, i.e., for large values of p, can be calculated by 
the free volume theory. 

The perturbation energy < U > HS (p) = Y.t=o h ^ {M v (x) - NM V (B%)) in Eq. 
(^D is given by the average values M v (x) of the Minkowski functional in the hard- 
sphere reference system. In general, the average values of the Minkowski functionals 
for correlated disks of radius R are given by 

M {x) = 1 - e-'X*' ™) , 

M l (x)= xu(x,cW)e- xf ^ c(n) \ (7) 
M 2 (x) = x (e{x, cW) - x (u(x, c< n ')) 2 ) e - xf ^ cM ) , 

where the functions f(x, c^), u(x, c^) and e(x, c^) depend on the normalized den- 
sity of disks x = pnR 2 and the hierarchy of correlation functions c^(xi, . . . , x n ) 



of the centers of the disks. In lowest order, i.e., for an approximation by a GauB- 
Poisson process one obtains the correlated average of the inter sectional area, bound- 
ary length, and integral curvature in terms of the two-point correlation function g(f) 
with c^(f) = (g(r) — l)p 2 : 



f(') = 1 - i Jdrc^ir) ( 1 - f arcsing _ iM^i _ 
u(x)= 1 - f / dfc^ (r) ( 1 - f arcsin (8) 
e(x) = l-±jdrcW{r)e{2R-\r\). 

P R 

For a Poisson distribution with g(r) = 1 one obtains f^(x,c^) = u^ p \x,c^ n ') = 
e( p >(x,c( n ') = 1. For a hard-core process approximated by the correlation function 
the mean values 

f \f\ < D 

9& = (9) 
I 1 \r\>D 

one obtains with 5 = D / (2R) < 1 



f^(R, D)= 1 + |(1- l5y/T=P{l + 25 2 ) - |(1 - 45 2 ) arccostf) 



u ( - hc \R,D)= l + x(l- l5^l^¥- f(l-2<5 2 )arccos5) 



>S hc \R,D)= l + 2x8 2 



From Eqs. (0), flf]), (0), and ([H]) it is possible to derive all thermodynamic properties 
of the system needed to construct the phase diagram. In the limit 5 = one recovers 
the equation of state for an ideal gas of overlapping disks. 



4 Monte- Carlo algorithm for overlapping disks 

With the definition @ of the Hamiltonian 7i(A) for a configuration A = Uj-B^(xj) 
of partially penetrable disks B^(xi) it is possible to perform directly Monte-Carlo 



simulations. Although lattice models based on the morphological Hamiltonian (Q) 
were defined and extensive computer simulations have been performed for complex 
fluids such as colloidal dispersions and microemulsions fl5f , not many continuum 
models are studied yet. A model based on correlated partially overlapping grains 
(see Figure P seems to be a promising starting point to study geometric features of 
complex fluids. In this section we describe the implementation of an algorithm in two 
dimensions for the parallel machine CM5 using the SIMD (single instruction multiple 
data) technique. 

In a first step, the Minkowski measures M U (A) are calculate analytically for a 
given configuration = Uf =l B^{xi) of N partially overlapping disks. Then, discs 
are added, removed, or shifted which leads to a new configuration A'. The difference 
in the Minkowski measures, i.e., in the Hamiltonian TC determines the probability 
for the acceptance of the new configuration A' according to the usual Metropolis 
dynamics. Unfortunately, the computational cost to evaluate the energy @ of a 
configuration given by Eq. (Q) is enormous so that an efficient algorithms is necessary 
in order to make Monte-Carlo simulations feasible. 

The Minkowski functionals M V (A N ) for the union {jf =1 B^{xi) may be calculated 
straightforwardly via the additivity relation (|^), i.e., as a sum of multiple overlaps 

M U (A N ) = Mv(Bi) — M u (Bi fl Bj) + . . . + (—l) N+l M v (Bx fl . . . fl B N ) . (11) 

i i<j 

which follows from Eq. (^) by induction. The right hand side of Eq. (|TT|) only 
involves convex sets and may be applied together with the definition of M v to compute 
M U (A N ). However, this algorithm becomes inefficient when the amount of overlap 
between the augmented balls is excessive, since one has to compute many redundant 
and mutually cancelling terms. Thus, this approach works only for 8 ~ 1. Therefore, 
we proceed alternatively as described in Ref. [[!(]]. The morphology of a configuration 
A (see Figure D) is unambiguous determined through the borderline between black 
and white regions. So it is possible to determine all Minkowski measures M v with an 



appropriate parameterization together with the local curvature of those borderlines. 
For instance, the two-dimensional area M can be calculated by applying Gauss' 
theorem 

M (A) = - [x-ndS (12) 



2 
OA 

where n denotes the normal vector to the boundary dA at x. Accordingly, one 
obtains the length of the boundary line M\{A) = J dS and the Euler characteristic 

dA 

M 2 (A) = Mi (A) I (7T.R) + Xi/(27r) as integrals along the boundary dA of a wetting 
layer shown in Figure [l], for instance. Here, \% denotes the angle between the two 
normals at the intersection point i of two disk boundaries. 

In usual Metropolis Monte-Carlo Simulations the initial state converges with a 
characteristic time scale r towards a configuration sampled from the equilibrium 
distribution. Whether or not such an equilibrium configuration is reached cannot be 
decided a priori and test runs have to be performed, in particular, for this model 
where the evaluation of the Boltzmann weights are difficult and extremely CPU-time 
expensive. In Figure || we show a typical time series N(t) of the number density of 
disks. At each sweep every disk were tried to move exactly once. During the first 
500 sweeps the systems relaxes from the initial configuration of 4800 disks to the 
thermal equilibrium of N = 4272 disks on average. The relaxation time is found to 
be t = 180 sweeps. At larger densities N the relaxation time increases considerably, 
so that equilibrium configurations could not be reached for densities x = pnR 2 > 3. 
Therefore, it was not possible to verify the existence of a second fluid- fluid phase 
transition at larger densities and a triple line predicted by a mean-field approximation 
of this model [14" . 

First, we use the algorithm to compute the average values of the Minkowski mea- 
sures M v (p) for (3 = as function of the density p of the particles. In Figure £|(a) the 
mean values of the covered area Mq (full circles), the length M\ of the liquid boundary 
lines (triangles) and the Euler characteristic M 2 (stars) are shown for three different 



ratios, namely 5 = (dashed lines), 5 = 0.4 (solid lines), and 5 = 1.0 (dot-dashed 
lines). Whereas for the pure hard-disk system (5 = 1.0) and for completely penetra- 
ble disks (5 = 0.0) a difference between the numerical results and theoretical values 
could not be found, one can observe deviations from the analytic expression (|TUp for 
large densities pmR 2 > 2 and 5 = 0.4. For small densities pnR 2 , the configuration 
A N consists of isolated disks with negligible overlap of the wetting layers; therefore, 
each measure starts out linearly in the number of particles. As the overlap increases, 
the volume Mq of the wetting layer saturates and the total boundary length Mi and 
M 2 of the coverage decreases. M 2 changes even its sign, since the curvature of the 
singular corners (vertices) at intersections of the wetting layers is negative and starts 
to dominate the positive contribution from the spherical parts of the boundary as the 
overlap increases (see Figure [l]). Negative values of the Euler characteristic M 2 are 
typical for a highly connected structure with many holes, i.e., undisturbed cavities of 
binary fluid. 

In Figure |3] the morphological measures of the wetting layer are shown as function 
of the chemical potential /i at finite temperature flhonR 2 = 1/0.355. The disconti- 
nuities in the measures indicate clearly a phase transition from a low density fluid 
phase (small volume M of the wetting layer) to a high density phase (large M ) 
where the colloidal particles move freely in a large region formed by the union of the 
wetting layers. In Figure [I] a configurational shnapshot close to the critical point at 
T c = l/(f3nR 2 h ) = 0.368 and p c 7rR 2 = 0.82 is shown, which illustrates both type of 
fluid phases. Note, that the typical size of the connected wetting layer (white region) 
is much larger than the distance between colloidal particles. 

The coexistence of the low-density fluid phase of isolated colloidal particles sur- 
rounded by a spherical wetting layer and the high-density fluid phase where the 
colloidal particles move freely in a fluid wetting phase is shown in Figure |](b). For 
the existence of the phase separation the inherent multi-particle interactions of the 
colloidal particles are crucial. For 5 = 0, for instance, a pair-potential approximation 



of the Hamiltonian (^j would not be bounded from below yielding unstable thermo- 
dynamic configurations. And for 5 << 1 only the multi-particle character of the 
Hamiltonian guarantees that for large densities p the interaction energy is limited so 
that the colloidal particles can move freely within the fluid film. 

Although the theoretical expression ([H]) for the Minkowski measures works well 
for densities x < 2, the phase coexistence calculated by perturbation theory is not in 
satisfying agreement with the simulation results. Since the main goal of this approach 
is an effective theory expressing the thermodynamic properties of the colloidal systems 
in terms of the thickness L(T, c) of the wetting layer, one needs to improve the 
perturbation theory before the analytic results (|10D can be applied to more realistic 
situations such as the experiments of silica beads in lutidine- water mixture [TJ, 0] . 

The evaluation of the Boltzmann factor is very CPU-time consuming which makes 
the simulation difficult, in particular, for large densities. An improved algorithm could 
verify the existence of a second fluid-fluid phase transition at larger densities which 
is indicated by perturbation theory due to multiple particle interactions of the Euler 
characteristic in the Hamiltonian. 
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Figure captions 



Fig.l: An ensemble of hard colloidal particles (black disks) surrounded by a 
fluid wetting layer (white ring). The interaction between the colloids is determined 
by the free energy of the fluid film (white region) which may cause a fluid-fluid 
phase separation of the hard particles. Thus, the spatial structure of the phases, 
i.e., the morphology of the white regions determines the configurational energy which 
determines itself the spatial structure due to the Boltzmann factor in the partition 
function of a canonical ensemble. A main feature of this model is the occurrence of 
different length scales: the clusters of the particles, i.e., the connected white regions 
are much larger than the 'microscopic' radius of the disks and the typical nearest- 
neighbor distance within a cluster. 

Fig. 2: Times series N(t) and corresponding correlation function c(t) where the 
first 3000 sweeps were omitted. The temperature T = 0.36 is slightly below the 
critical point at T c = 0.368 for h — 1, h± — h 2 — 0, where T = kBT/(nR 2 h ) denotes 
the reduced temperature. 

Fig. 3: Isotherms at f = 0.355 {5 = 0.4, h = 1, hi = h 2 = 0) of the density 
p, the covered area F = M , the length U = 2nMi of the liquid boundary lines and 
the Euler characteristic x — measuring the topology of the configurations. One 
can clearly observe the phase transition at @fi = 2.1. The functional form of the 
morphological measures are typical: x exhibits a discontinuous jump from positive 
to negative values, U has its maximum at the transition, whereas F is a monotonous 
increasing function. 



Fig. 4; (a) Isotherms at (3 = 0.0 and 5 = 0.4 of the covered area Mo (full circles), 
the length Mi (triangles) of the liquid boundary lines and the Euler characteristic M 2 
(stars). The analytic result for 5 = 0.4 (solid lines) given by Eq. (0) works quite well 
for densities x < 2. For 5 = 0.0 (fully penetrable disks, dashed lines) and for 5 = 1.0 
(hard disks, dot-dashed line) no differences between numerical determined values 
and theoretical predictions could be found, (b) Phase diagram with hi — hi — 0, 
h = 1, and 5 = 0.4; m = tiR 2 . The dashed line indicates the location of the 
coexisting densities. The perturbation result (solid line) based on the isotherms (a) 
overestimates the stability of the homogeneous fluid phase yielding a lower critical 
temperature. 
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